%%
%clear; close all;clc;
%%
vardata = readtable('allVarData.xlsx','Sheet',char(clabel(country)));
jln_ui = readtable('IndexQuarterlyData.xlsx','Sheet',char(clabel(country)));
%cedata = readtable('CE gdp fcasts update.xlsx','sheet','fcasts');
bbddata = readtable('BBD indices update.xlsx','sheet','bbdvar');
%globaldata = readtable('IndexQuarterlyData.xlsx','Sheet','Global'); 

dates=table2array(jln_ui(:,1));
datevector =datevec(dates);
jlndates=datestr(datenum(datevector(:,1),datevector(:,2),1));
jln_ui= [cellstr(jlndates) jln_ui(:,2:end)];
jln_ui.Properties.VariableNames{'Var1'} = 'date';

%bbd uncertainty data
bbddates = table2array(bbddata(:,1));
datevector = datevec(bbddates);
bbddates=datestr(datenum(datevector(:,1),datevector(:,2),1));
bbddata(:,1)=cellstr(bbddates);
bbddata = bbddata(:,[1 country+1]);
bbddata.Properties.VariableNames{char(clabel(country))} = strcat(char(clabel(country)),'_bbd');
bbddata.Properties.VariableNames{'Var1'} = 'date';

%Var data
vardates=char(table2array(vardata(:,1)));
keySet = {'Q1','Q2','Q3','Q4'};
valueSet = [3 6 9 12];
datedict = containers.Map(keySet,valueSet);
date_part=char(vardates);
quarter=date_part(:,1:2);
year=str2num(date_part(:,4:end));
varDates=[];
for t=1:length(vardates)
    nextDate = datestr(datenum(year(t),datedict(quarter(t,:)),1));
    varDates = [varDates; nextDate];
end
vardata(:,1) = cellstr(varDates);
varlength=1:size(vardata,1);
varlength=varlength';
vardata(:,9) = array2table(varlength);
vardata.Properties.VariableNames{'Var9'} = 'Index';
if ismember('Var1',vardata.Properties.VariableNames)
    vardata.Properties.VariableNames{'Var1'} = 'date';
end

%1.merge in jln indices
testdata = outerjoin(vardata,jln_ui,'MergeKeys',true);
testdata=sortrows(testdata,'Index');

if muindex==0
    %3.Merge in BBD index
    testdata = outerjoin(testdata,bbddata,'MergeKeys',true);
    testdata=sortrows(testdata,'Index');
end

dates=datenum(table2array(testdata(:,1)));
data=table2array(testdata(:,2:end));
varnames = testdata.Properties.VariableNames;
varnames(:,1)=[]; % drop date
data(:,8)=[];varnames(:,8)=[];  %drop index used in join

%Drop NaNs
missing_values = any(ismissing(data),2);
data=data(~missing_values,:);
dates = dates(~missing_values);
disp([datestr(dates(1)) ' ' datestr(dates(end))])

clearvars -except data varnames dates projectdir projectdirDB clabel country NumberConsecutiveshocks LogVsLevel muindex ...
                  nn ydk zdk


if LogVsLevel
    if muindex==0
logEPU_rescaled =     data(:,find(strcmp(varnames,strcat(clabel{country},'_bbd')))); %it is already logged in Chrsi Redl data
    else
logEPU_rescaled =     data(:,find(strcmp(varnames,'Fin_Index')));
    end
else
    if muindex==0
logEPU_rescaled = exp(data(:,find(strcmp(varnames,strcat(clabel{country},'_bbd')))));
    else
logEPU_rescaled = exp(data(:,find(strcmp(varnames,'Fin_Index'))));        
    end
end
% pos = [find(strcmp(varnames,'Credit')) ...
%        find(strcmp(varnames,'BankRate')) ...
%        find(strcmp(varnames,'Hours')) ...
%        find(strcmp(varnames,'GDP'))];
if sum(ismember(varnames,'Hours'))
pos = [find(strcmp(varnames,'Credit')) ...
       find(strcmp(varnames,'BankRate')) ...
       find(strcmp(varnames,'Hours')) ...
       find(strcmp(varnames,'GDP')) ...
       find(strcmp(varnames,'Investment')) ...
       find(strcmp(varnames,'CPI'))];   
else
if country ~= 7 
pos = [find(strcmp(varnames,'Credit')) ...
       find(strcmp(varnames,'BankRate')) ...
       find(strcmp(varnames,'Employment')) ...
       find(strcmp(varnames,'GDP')) ...
       find(strcmp(varnames,'Investment')) ...
       find(strcmp(varnames,'CPI'))];       
else
pos = [find(strcmp(varnames,'Credit')) ...
       find(strcmp(varnames,'DiscountRate')) ...
       find(strcmp(varnames,'Employment')) ...
       find(strcmp(varnames,'GDP')) ...
       find(strcmp(varnames,'Investment')) ...
       find(strcmp(varnames,'CPI'))];           
end
end
%{'CPI'}  {'Hours'} {'Investment'}  {'Consumption'}  {'GDP'}  {'Credit'} is in log levels 
x_ = data(:,pos);
 
 dataMatrix= 100*[logEPU_rescaled/100               x_(:,2)/100 log(x_(:,4)) log(x_(:,3))];
 dataName   = {'EPU',               'Short-term nominal rate','Output','Employment'}; 
%dataMatrix= 100*[logEPU_rescaled/100  log(x_(:,1)) x_(:,2)/100 log(x_(:,4)) log(x_(:,3)) log(x_(:,5:6))]; 
%dataName   = {'EPU','Credit index','Short-term nominal rate','Output','Employment',  'Investment','CPI'}; 

clear ss_ se_ x_

% 

pos_shock = 1;

nobs       = size(dataMatrix,1);

%-------------------------------------------------------------
% Step 0: parameter settings
%-------------------------------------------------------------
nlags  = 2;                   % number of lags in VAR  
nconst = 1;                   % 3 = for constant, linear and quadratic trend, 2 = for constant and trend, 1 = for constant, 0 otherwise 

%-------------------------------------------------------------------------%
prior=1;                      % for prior,  1   = standard prior
                              %             inf =     flat prior

%set prior for first lag of each variable (if estimation with Minnesota prior)
[~,coly]=size(dataMatrix);
lev=ones(coly,1);                                          
%-------------------------------------------------------------------------%       
  nimp = 36;                  % horizon of VD and IRFs to be shown         

est_method = 0;               % 0 : OLS point estimation
                              % 1 : Bayesian estimation with Minnesota prior 
                              % 2 : OLS point estimation and classical bootstrap    
                              
%-----------------------------------------------------------------------
% Step 1: load and transform data
%-----------------------------------------------------------------------
[y,x] = mk_vdm(dataMatrix,nlags,0);
%y is T x nvars matrix of rhs variables
%x is T x nvars*nlags  of lhs variables
[T,nvars]=size(y);

%---------------------------------------------------------------------
% Step 3: Estimate VAR, extract shocks, FEV shares explained, and IRFs 
%---------------------------------------------------------------------
if nconst==1
   const = ones(T,1); 
   x = [x const];
   ncoeffs = nvars*nlags+1;
elseif nconst==2
   x = [x ones(T,1) (1:T)'];
   ncoeffs = nvars*nlags+2;
elseif nconst==3
    x = [x ones(T,1) (1:T)' ((1:T).^2)'];
%   dummy_1975Q2 = dates_q_==1975.5;
%   x = [x ones(T,1) (1:T)' ((1:T).^2)' dummy_1975Q2(1:end-nlags)];
   ncoeffs = nvars*nlags+3;   
else
   ncoeffs = nvars*nlags; 
end

%ncoeffs + 1 x nvars matrix of regression coefficients
b=x\y;       %OLS coefficients: ncoeffs x nvars (** more efficient way of computing least squares; i.e. b=inv(x'*x)*x'*y;
res = y - x*b;          %T x nvar matrix of residuals
vmat = (1/T)*(res'*res);
stds=sqrt(diag(vmat));

%OLS point estimates
[FEVD_IRF,v]  =RecursivenessIdentification(b,res,vmat,nvars,nlags,nimp, pos_shock);

%----------------------------------------------------------------------
% Step 3: Plots
%----------------------------------------------------------------------
% clearvars -except dataName VarsToPlot FEVD_IRF plottype nvars
% save(strcat('IRFDATA_',plottype))

% %FEV shares explained by shock
% plotFEV(FEVD_IRF,vars,slope,use_yields,shock_extract);
% 
%IRFs to shock
% plotIRF(FEVD_IRF,vars,slope,use_yields,shock_extract)
    %plot(FEVD_IRF(:,[nvars+posGDP])*100)
uEPUVAR = [NaN(nlags,1); v];
clearvars -except NumberConsecutiveshocks logEPU_rescaled uEPUVAR FromM2Q dataMatrix nlags dates country muindex clabel ...
           nn ydk zdk 


%% FIT AR(p) as in Brogaard and Detzel (MS, 2015)

%% FIT AR(1)  --> Shocks 1961/M1 -- 2019/M12
%dates_EPU = dates_EPU(13:end);

%%
    uEPU = uEPUVAR  ; clear  uEPUAR1 uEPUAR4            uEPUVAR 
%   NumberConsecutiveshocks = 2;
% % NumberConsecutiveshocks = 4;
    
    uEPU = uEPU/nanstd(uEPU);
clearvars -except NumberConsecutiveshocks logEPU_rescaled uEPU dataMatrix nlags dates volaShock country muindex clabel ...
                  nn ydk zdk
% /////////////////////////////////////////////////////////////////////////
%----------- Construct indicators for seq. of shocks with 2 ident schemes-%
IndPreviousNShocksEPU= zeros(size(uEPU));
% /////////////////////////////////////////////////////////////////////////
for tt=NumberConsecutiveshocks:length(uEPU)-1 
    if tt==NumberConsecutiveshocks
%      %if (                                      sum(uEPU(tt-NumberConsecutiveshocks+1:tt)>0)==NumberConsecutiveshocks && uEPU(tt+1)<0);  IndPreviousNShocksEPU(tt) = 1; end
        if (                                      sum(uEPU(tt-NumberConsecutiveshocks+1:tt)>0)==NumberConsecutiveshocks                );  IndPreviousNShocksEPU(tt) = 1; end
    else
%      %if (uEPU(tt-NumberConsecutiveshocks)<0 && sum(uEPU(tt-NumberConsecutiveshocks+1:tt)>0)==NumberConsecutiveshocks && uEPU(tt+1)<0);  IndPreviousNShocksEPU(tt) = 1; end
%       if (                                      sum(uEPU(tt-NumberConsecutiveshocks+1:tt)>0)==NumberConsecutiveshocks                );  IndPreviousNShocksEPU(tt) = 1; end
        if (uEPU(tt-NumberConsecutiveshocks)<0 && sum(uEPU(tt-NumberConsecutiveshocks+1:tt)>0)==NumberConsecutiveshocks                );  IndPreviousNShocksEPU(tt) = 1; end
    end
end
clear tt